Ab initio study of gap opening and screening effects in gated bilayer graphene 
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The electronic properties of doped bilayer graphene in presence of bottom and top gates have 
been studied and characterized by means of density-functional theory (DFT) calculations. Varying 
independently the bottom and top gates it is possible to control separately the total doping charge 
on the sample and the average external electric field acting on the bilayer. We show that, at fixed 
doping level, the band-gap at the K point in the Brillouin zone depends linearly on the average 
electric field, whereas the corresponding proportionality coefficient has a nonmonotonic dependence 
on doping. We find that the DFT-calculated band-gap at K, for small doping levels, is roughly 
half of the band-gap obtained with standard tight tinding (TB) approach. We show that this 
discrepancy arises from an underestimate, in the TB model, of the screening of the system to the 
external electric field. In particular, on the basis of our DFT results we observe that, when bilayer 
graphene is in presence of an external electric field, both interlayer and intralayer screenings occur. 
Only the interlayer screening is included in TB calculations, while both screenings are fundamental 
for the description of the band-gap opening. We finally provide a general scheme to obtain the full 
band structure of gated bilayer graphene for an arbitrary value of the external electric field and of 
doping. 



PACS numbers: 71.15.Mb, 73.22.-f, 73.61.-r, 81.05. Uw 
I. INTRODUCTION 

Among the nanoscale forms of carbon, bi- 
layer graphene has recently attracted much 
interest P2 B | 4 | 5 | 6 | 7 | 8 | Indeed, it has been found, both 

theoretically and experimentally, that in presence of an 
asymmetry between the two graphene layers, generated 
by an external electric field, a band-gap can be opened. 
This makes bilayer graphene a tunable-gap semicon- 
ductor and therefore an exciting structure for future 
application in nanoelectronics. 

In particular, in the experiments of Ohta et alP bi- 
layer graphene is synthesized on silicon carbide (SiC) 
substrate, and a small n-type doping is acquired by the 
system from the substrate. In this case, the bilayer sym- 
metry is broken by the dipole field created by the deple- 
tion of charge on SiC and accumulation of charge on the 
bilayer. Additional n doping is induced by deposition 
of potassium atoms above the bilayer. Varying the con- 
centration of potassium atoms, one can vary the asym- 
metry between the two sides of the system and measure 
the electronic properties and the band-gap opening by 
angle- resolved photoemission spectroscopy (ARPES). 

Oostinga et alP used a double-gated system, where 
monolayer and bilayer graphenes are placed between two 
dielectrics, which act as bottom and top gates. The 
double-gated structure gives the possibility to control in- 
dependently the doping level and the perpendicular elec- 
tric field acting on the system. In this configuration, 
they measure the dependence of the resistance on the 
temperature and on the electric field. They observe a 
gate-induced insulating state in bilayer graphene which 
originates from the band-gap opening between the va- 
lence and conduction bands. 



As for theoretical studies, McCamPused a tight bind- 
ing (TB) model to study the band structure of the bilayer 
graphene in presence of an energy difference between the 
two layers, which determines a band-gap opening. In 
particular, he considered a single gate acting on the sys- 
tem, and he found a roughly linear relation of the gap 
with the accumulated charge n on the bilayer, for n val- 
ues up to 10 x 10 12 cm~ 2 . Min et alP performed ab initio 
density-functional theory (DFT) calculations of undoped 
bilayer graphene in a constant external electric field, us- 
ing the generalized gradient approximation (GGA) for 
the exchange-correlation functional. They confirmed the 
general picture provided by the TB model, although DFT 
screening results stronger. Moreover, Aoki et AmawashP 
performed an ab initio DFT study on the band structure 
dependence of undoped layered graphene on the stacking 
and external field, using the local density approximation 
(LDA) for the exchange-correlation functional. In con- 
trast with the GGA study of Min et al.pl their results on 
undoped bilayer graphene in a uniform external electric 
field are in agreement with the TB ones. Again, Castro 
et al.EEl showed experimentally, by measuring the Hall 
conductivity and the cyclotron mass of biased bilayer 
graphene, and theoretically, using TB methods, that a 
band-gap between the valence and conduction bands can 
be tuned by an applied electric field. 

Other theoretical DFT studies have been devoted to 
the understanding of the band structure dependence on 
the stacking geometry,^ on the presence of adsorbed 
molecules p2 and to the analysis of the distribution of 
the induced charge densities.^ Instead, other experi- 
mental st udies focu sed on the Raman spectra of bilayer 
graphene.^ 14 ' 15 1 16 1 Recently, an experimental work on in- 
frared spectra of gated bilayer graphene as a function of 
doping appeared,^ and a comparison with TB calcula- 
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tions suggests that the TB prediction of the gate-induced 
band-gap is overestimated. 

In this work, we study by means of DFT ab initio cal- 
culations the band-gap opening in bilayer graphcne, both 
as a function of the external electric field and as a func- 
tion of doping. The paper is organized as follows: in Sec. 
[IT] a description of the system we investigate is reported, 
along with the computational details. Results are pre- 
sented in Sec. 
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where the dependence of the gap on 
the electric field and doping is first shown and compared 
with TB calculations. Then, a detailed analysis of the 
screening properties of the bilayer to the external electric 
field is reported. The effect of the electronic temperature 
on the screening is also investigated, and the nonmono- 
tonic behavior of the band-gap as a function of doping 
at fixed electric field is explained. The GW correction of 
the DFT-calculated response of the gap to the external 
electric field is presented. Finally, a general scheme to 
obtain the full band structure of gated bilayer graphcne 
is provided, and a comparison with experimental findings 
is reported. In Sec. |IV| our conclusions are drawn. In the 
Appendix we describe in detail the top and bottom gates 
implementation in our DFT calculations. 



II. 



THEORETICAL BACKGROUND 



A. Bilayer graphene in bottom and top gates 

The experimental setup where bilayer graphene feels 
different bottom and top gates is schematically repre- 
sented in Fig[T] The bilayer is first grown on a dielectric 
material, of width D2 and relative dielectric constant t r i- 
Applying a voltage difference V g 2 (bottom gate) between 
the dielectric and the bilayer, a doping charge per unit 
surface nie — £2^2 is accumulated on the bilayer, where 
e is the electron charge (e = - |e|) and £2 = ^r2/D 2 - 
eo is the permittivity of the vacuum. Depositing an- 
other dielectric material of width Di and with relative 
dielectric constant e r \ over the bilayer, and applying a 
gate voltage V g \ (top gate) between the dielectric and 
the bilayer, an additional doping charge per unit sur- 
face ri\e = £,iVgi is accumulated, where £1 = e a e r i/Di. 
A total doping charge ne is therefore accumulated on 
the bilayer, where n = n\ + According to standard 
notation, positive n corresponds to electron doping and 
negative n corresponds to hole doping, p\ and p2 are 
the electronic charges per unit area (with respect to the 
neutral case) accumulated on layer 1 and layer 2, respec- 
tively. In particular, the sum of p\ and P2 is determined 
by the electrostatics, and it is equal to the sum of n\ and 
7i2- However, the individual values of p\ and p 2 depend 
on the screening properties of the system, and in general 
Pi 7^ ni and p 2 7^ n 2 . In the configuration shown in 
Figjl] layer 1 and layer 2 of the bilayer feel an electric 
field Ei and E2, respectively, given by 
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FIG. 1: Schematic representation of the experimental setup 
where bilayer graphene is placed between two dielectric ma- 
terials, and it is doped by applying bottom (V g 2) and top 
(V^i) gates. The width of the two dielectrics is much larger 
than the distance between the dielectrics and the bilayer. A 
doping charge per unit area ne — (m + n 2 )e is accumulated 
on the bilayer, where ii2e = £ 2 V^ 2 and n\e = £iV 3 i are the 
charges from bottom and top gates, respectively, pie and p 2 e 
are the electronic charges per unit area (with respect to the 
neutral case) on layer 1 and layer 2, respectively. In partic- 
ular, (pi + p 2 )e = (m + n 2 )e, while pi / n\ and p 2 / n 2 . 
Layer 1 and layer 2 of bilayer feel electric fields E\ = —nie/eo 
and E2 = n 2 e/eo, which determine an average electric field 
-Bav = (ri2 — ni)e/(2eo). eo is the permittivity of the vacuum. 



E 2 = n 2 e/e . 
The average electric field E av is defined as 

E av = (E 1 + E 2 )/2 

= (n 2 - ni) e/(2e ) 
= K-n 2 ) |e|/(2e„). 



(2) 



(3) 



£1 



-ni e/e , 



(1) 



Positive E av is oriented from dielectric 1 to dielectric 2 
(i.e., from top to bottom gate). When n\ and n 2 are 
equal, we are in the case of equal bottom and top gates, 
and E av vanishes. When n\ is zero, we are in presence 
of bottom gate alone. The top gate can also be gener- 
ated by a chemical doping, with the deposition of alkali 
or halogen atoms on the bilayer. In this work, the elec- 
tric fields E\ and E 2 are simulated using periodically 
repeated boundary conditions by introducing dipolc and 
monopole potentials, as described in the Appendix. 

The presence of different bottom and top gates gener- 
ates an electrostatic potential which is different on layer 
1 with respect to layer 2, and this asymmetry gives origin 
to a band-gap opening. In Fig(2]we show the band struc- 
ture of undoped bilayer graphene, in absence of bottom 
and top gates, where no gap is observed, and in presence 
of different bottom and top gates in which case a gap is 
opened. In order to simplify the discussion, in this work 
we define a signed gap U at the K point in the Brillouin 
zone (BZ), which is negative for E av < 0, and positive 
for E w > 0. 
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FIG. 2: Band structure, around the K point in the BZ, of 
undoped bilayer graphene in absence of bottom and top gates 
(solid line) and in presence of different bottom and top gates 
(dashed line). 



FIG. 3: Schematic representation of the electrostatic model 
used with TB calculations. The bilayer of thickness d feels an 
external electric field E\ on layer 1 and E2 on layer 2, which 
induce a charge per unit surface pie and p2e on layer 1 and 
layer 2 and a local electric field E* . 



B. Computational details 

The presented ab initio results based on the DFT, 
are done using both the Perdew-Burke-Ernzerhof (PBE) 
(RefP) GGA, and the Perdew-Zunger (PZ) Ref.^ 
LDA exchange-correlation functionals. Core elec- 
trons are taken into account using the pseudopo- 
tential method, with norm-conserving Troullier-Martin 
pseudopotentials.^ Plane-waves basis set is used to de- 
scribe valence electron wave functions and density, up 
to a kinetic energy cut-off of 40 and 600 Ry, respec- 
tively. The electronic eigenstates have been occupied 
with a Fermi-Dirac distribution, using an electronic tem- 
perature of 300 and 30 K. The BZ integration has been 
performed with a uniform k point grid of (80 x 80 x 1) 
and (240 x 240 x 1) for the two temperatures, respec- 
tively. The experimental lattice constant a = 2.46 A of 
two-dimensional graphite is used. The layer-layer dis- 
tance d is fixed at the value of 3.35 A, as in graphite. 
The two layers are arranged according to Bernal stack- 
ing. The length L of the supercell along z is 17.2 A. Cal- 
culations have been performed usine ; the PWscf codecs of 
the Quantum ESPRESSO distribution^. 

In this work we perform also TB calculations. The TB 
model we use is characterized by two parameters, 7y and 
7j_, which represent the first nearest-neighbors in- plane 
hopping and the interplane hopping between vertically 
superposed atoms in the Bernal stacking configuration, 
respectively. 7y is related to the Fermi velocity in sin- 
gle layer graphene, vj — 7||av3/(2?i). We use a value 
of 711= 3.1 eV, as inferred from experimental measure- 
ment IM3I23I Within the TB model, 27^ corresponds to 
the band splitting between the lowest occupied tt band 
and the highest unoccupied tt band at K (see Fig{2]). We 
use a value of 7j_ = 0.4 eV , as obtained from experi- 
mental measurements G0I211ZI Similar values for these TB 



parameters have been used in literature! 5 * 7 * 8 ! In addition 
to these parameters, we consider the energy difference 
between layer 2 and layer 1 induced by the electric field, 
which coincides with the signed gap U at the K point in 
the BZ (see Fig{2]and RefP). 

In the TB formalism, in order to obtain the relation 
between the gap U and the average electric field E av , 
a simple electrostatic model is used. In Fig(3] we show 
a schematic representation of this electrostatic model. 
Charges per unit surfaces p\e and p 2 e are concentrated 
on the two layers of the bilayer, which create a screened 
field E* inside the system. Using simple electrostatic 
equations, we have 



E* = E, 



Ap e 



(4) 



where E av = {E\ + E 2 )/2 and Ap = p 2 — pi- Ap is 
calculated from the square modulus of the eigenfunctions 
in the two layers. The energy difference between layer 2 
and layer 1, i.e., the band-gap U at K, is given by 



U = - d E* 



(5) 



Inserting into eq.([5]) the expression of E* as given in 
eq.Q, and writing Ex and E 2 as in Eq.Q and ([2]), we 
obtain: 



U 



de? 
2e 



(ni - n 2 + Ap). 



(6) 



Therefore, in the TB calculations the electronic screening 
is evaluated using the simplified electrostatic model de- 
scribed above, contrary to the DFT formalism where the 
detailed shape of the charge distribution is fully taken 
into account. 
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FIG. 5: (Color online) a as a function of doping n, for an 
electronic temperature of 300 K, calculated with DFT-GGA 
(x crosses) and LDA (+ crosses), calculated using TB model 
with 7|| = 3.1 eV and 7x = 0.4 eV (circles), and using GW 
correction (up-triangles). The GW correction is obtained as 
described in Sec. |III E| The continuous thick (black) line is 
the fit, as written in Eq.(33l, of the GW result. The value 



of a in absence of electronic screening is independent on the 
doping, and it is a barc = de 2 /(2e ) = 30.3 x 10~ 12 cm 2 meV. 



experimentally obtained in bilayer graphene by the appli- 
cation of a gate voltage with a SiC>2 dielectric^ or with 
a polymeric electrolyte.^ Our results show that U has a 
linear dependence on the applied electric field E av . We 
therefore define a linear response a(n) such that 



FIG. 4: (Color online) DFT-GGA calculated U as a function 
of (711—712), i.e., the average electric field divided by |e|/(2eo), 
for electron and hole dopings. Two values of electron doping 
are shown, n = 38 and 2 x 10 12 cm -2 , and two values of hole 
doping, n — -38 and -2 x 10 12 cm -2 . Points are the calculated 
values, for an electronic temperature of 300 K, while lines are 
linear fits. 



U (n,E av ) = a(n)(n x - n 2 ). 



(7) 



III. 



RESULTS 



Band gap as a function of the external electric 
field and doping charge 



As anticipated in Sec jll A[ when bilayer graphene feels 
different bottom and top gates, a band-gap U is opened. 
In this section we first investigate the dependence of U 
on the average external electric field E av , at fixed doping 
n. 

In Figj4] we show the DFT-GGA calculated U as a 
function of [n\ — 712) [i.e., the average electric field E av 
divided by |e| / (2e )] , for two values of electron and hole 
dopings. These values of doping are chosen as repre- 
sentative of two different doping regimes, which can be 



In Fig|5] we show a as a function of doping n, calcu- 
lated from Eq.(|7|, for an electronic temperature of 300 
K, within the DFT-GGA and LDA functionals, and using 



the TB model described in Sec II B Contrary to previ- 
ous results in literature our LDA and GGA results are 
very similar, and not in agreement with the TB ones. In 
particular, our results for zero doping and GGA func- 
tional are in agreement with the previous GGA study. 6 
They disagree instead with the ones computed with LDA 
functional in RefP. This is probably due to the fact 
that in RefP the authors used a coarse k point sam- 
pling (10x10x1) with respect to the ones used in this work 
and in RefP, and their results are likely unconverged. In 
the following we only present our GGA results. Both 
DFT and TB a's display a nonmonotonic behavior as 
a function of doping. However, the values of the DFT- 
calculated a are substantially different from the TB ones, 
especially for low doping values, i.e., when the Fermi level 
is close to the band-gap edges. This is the most interest- 
ing case for the application of the bilayer as active device 
in electronics. 

We notice that in the absence of electronic screening, 
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a is independent of the doping, and it is 



2e n 



= 30.3 x KT 12 cm 2 meV. 



(8) 



Thus, with the inclusion of the electronic screening, the 
DFT-calculated a becomes roughly three times smaller 
than the cv bare , which suggests that the screening effects 
are crucial for the description of the band-gap. 

In order to understand the origin of the difference be- 
tween the DFT and TB results, we notice that this can 
be due (i) to the calculated electronic band structure and 
charge transfer and (ii) to the electrostatic model used 
in the TB calculations, which gives a simplified descrip- 
tion of the crucial screening effects, fully included in the 
DFT formalism. To verify the quality of the electrostatic 
model, we introduce the quantity 77(71) defined as: 



Ap (n,U) = T)(n)U. 
Ap = pi — p\ is calculated from 



00 

i>2 = I [p{z) - po(z)] dz, 

'0 
M 



Pi 



[p(z) - po(z)} dz, 



(9) 

(10) 
(11) 



where p(z) is the planar average of the electronic charge 
density (per unit volume) for a doping n and in presence 
of E av and po{z) is the planar average of the electronic 
charge density (per unit volume) for the neutral case, 
with E av =0. Here and in the following z — indicates 
the plane at the midpoint of the two graphene layers. In 
our DFT calculations, ±00 corresponds to ±L/2, where 
L is the length of the supercell along z. 

77 is a measure of the charge transfer between layers in 
the presence of a band-gap U. We introduce this quantity 
because it is a direct outcome of the TB calculations, 
and no further electrostatic model is needed to compute 
it. Moreover, according to the electrostatic model used 
together with the TB formalism, described in Sec. II B 



the relation which gives a as a function of 77 is obtained 
dividing Eq.([6j| by U, 



oi(n) 



.hare 



1 - 77(71) a bare ' 



(12) 



In Fig(6] we show 77(71), calculated from Eq.(|9|, for an 
electronic temperature of 300 K, within the DFT and us- 
ing the TB model described in Sec |IIB| 77 as a function 
of n has a nonmonotonic behavior as found for a(n), and 
this trend is well described by the two methods. How- 
ever, the values of 77 calculated with the two formalisms 
arc different. Since no electrostatic model is used in the 
TB calculations, we conclude that this discrepancy orig- 
inates only from the difference between the DFT- and 
TB-calculated band structure and charge transfer. 

Moreover, comparing the DFT and TB results of a(n) 
and 77(77), we see that, for low doping levels, the relative 
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FIG. 6: (Color online) 77 as a function of doping n, calculated 
with DFT-GGA (crosses), with TB model (circles), and with 
GW correction (up-triangles), for an electronic temperature 
of 300 K. The GW correction is obtained as described in Sec. 

irrrEi 



difference, with respect to DFT values, of the TB/DFT 
a(n)'s is around 60%, while the analogous difference for 
77 is around 30%. Therefore, the electrostatic model used 
to compute a(n) in the TB formalism introduces a large 
error in the description of the band-gap opening in pres- 
ence of an external electric field. 



B. Electronic screening effects 

In this section we analyze where the simplified elec- 
trostatic model used in the TB calculations fails, and 
we propose a more sophisticated one. First of all, we 
know that in the TB formalism the energy difference be- 
tween the two layers coincides with the band-gap U at 
the K point in the BZ. In Fig(7]we show the band-gap 
U as a function of V2-1 — Vi ~ V\, where V% and V\ are 
the planar average of the DFT-calculated ionic, Hartree, 
and electrostatic potential energy on layer 2 and layer 1, 
respectively. The inclusion of the exchange-correlation 
potential does not change the result. We can notice that 
even in DFT formalism, U is correlated with the po- 
tential energy difference between the two layers, and in 
particular, 



(13) 



where (3 — 1.072, slightly higher than the expected uni- 
tary slope. 

To better understand the screening effects in the sys- 
tem, we investigate the linearly induced charge (per unit 
volume) p*- 1 -* 



pW (z;n,E av ) 



dp(z;n,E a 



E„ 



(14) 
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FIG. 7: DFT-GGA calculated U as a function of the differ- 
ence between the planar average of the ionic, Hartree, and 
electrostatic potential energy on layer 2 and layer 1, V2-1 , 
for different values of doping and of _B av . Points are the cal- 
culated values; the line is the linear fit, U = 1.072 x V2-1. 



viously no charge transfer between layers occurs, and the 
electronic screening to the external electric field is only 
characterized by an intralayer polarization. Moreover, 
we notice that the dependence of the induced charge on 
the doping is negligible. 

In Figl&ta) we show for bilayer graphene in pres- 
ence of the same external electric field. First of all, we 
notice that in the monolayer and in the bilayer are of 
the same order of magnitude. Then, we observe that the 
electronic screening of the bilayer to the external electric 
field, is characterized by (i) a charge transfer between 
the two layers, which is peculiar to the bilayer and (ii) 
an intralayer polarization, which is also present in the 
monolayer. 

In order to separate in the bilayer the interlayer from 
the intralayer polarization, we notice from Fig(8]that the 
intralayer induced charge is antisymmetric with respect 
to each individual layer. Thus we decompose the induced 
charge in the bilayer into a symmetric component, , 
and an antisymmetric component, pcP, with respect to 
each individual layer, pi 1 ^ and p^ are defined for z S 
{— d; d}, i.e., in an interval of width d around each layer, 
where d is the intralayer distance; they are calculated as 
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FIG. 8: Planar average of the linear induced charge (per 
unit volume) p"' of a graphene monolayer in presence of an 
external electric field i? a v 
doping level n — 1 x 10 12 c 
~ 2 (dashed line). 



1.6 x e/(2e ) 10 12 cm -2 for a 
2 (continuous line) and n = 19 



[p(z; n, E av ) - p(z; n, -E av )} 



(15) 



P?/.(*) = \ ±P (1) Mg n W d ~ *]} ■ (16) 



The symmetric, ps , and antisymmetric, Pa' , compo- 
nents are related to the charge transfer between the two 
layers and to the intralayer polarization, respectively. 

In Figjsjb) we show the symmetric component p^\ 
with respect to each layer, of the induced charge p^ 
shown in Fig(9ja). In Figj&J-c) we show the antisymmetric 
component pa . In particular, is very similar to the 
induced charge in the monolayer (Fig|8|, and it is of the 
same order of magnitude of the total induced charge in 
the bilayer [Fig{9 a] . On the basis of this qualitative anal- 
ysis of the linearly induced charge, we conclude that the 
intralayer polarization, which is not taken into account 
in the TB formalism, gives an important contribution to 
the screening properties of the system. 

In order to quantify the effect of the induced charge on 
the gap, we write the exact expression of the potential 
energy difference V2-1 in terms of the linearly induced 
charge and of the external average electric field E av 
using the Poisson equation in one dimension. We obtain 
the following: 



where p(z; n, -Eav) is the planar average of the charge den- 
sity (per unit volume) at a given doping level n and 
in presence of an external average electric field E av . 
Such p^ is antisymmetric with respect to z—0, i.e., 
p^ 1 \z;n,E ecv ) — — p*- 1 ^— z; n, E av ). In our plots we use 



the finite difference expression of p^\ i.e., Eq.(15l 



In Figj8]we show p' 1 ' for the graphene monolayer in 
presence of an external electric field E av — 1.6 x e/(2eo) 
10 12 cm -2 for two different doping levels. In this case, ob- 



V2-1 = V(d/2) - V(-d/2) 



de E a 



.2 r +oo d 



\^-z\p^(z)dz + 

2 p-\-oq j 

+ 2V J x |-2~ z|p(1)(z)dz ' (17) 

where ± d/2 = ± 1.675 A is the z coordinate of the two 
layers. Considering that p^(z) — —p^{—z), by simple 
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FIG. 9: a) Planar average of the linearly induced charge 
(per unit volume) p^ 1 ' for bilayer graphene in presence of an 
external electric field -E av = 1.6 x e/(2eo) 10 12 cm~ 2 for a 
doping level n — 2 x 10 12 cm -2 (continuous line) and n 
38 x 10 12 cm -2 (dashed line); b) symmetric component, p. 
and c) antisymmetric component, p„ , with respect to each 
layer, of the linearly induced charge p^ 1 ' shown in a) for the 
same doping levels. 



(i). 



algebra and without approximations we can rewrite V2-1 



FIG. 10: d a and d s as a function of doping n, as defined in 



Eqs.(22l and (231 



as 



V2-I 
where 



ed E a 



2e 



Ap + e D a ~ e D s , (18) 



D a 



\z~h tiKz) dz, 
£o Jo 1 



and 



e 

e o Jo 

Ap= I pM{z) dz - 



\ z --\pM(z) dz, 



p {1 \z) dz. 



(19) 
(20) 

(21) 



D a and D s represent the contributions to the poten- 
tial energy difference V2-1 given by the antisymmetric 
and symmetric components of the linearly induced charge 
around each layer. Equation (18) gives the exact expres- 



sion of V2-1 as a function of the external electric field 
and of the screening charge. 

We now rewrite D a and D s as follows: 



D a = d a {n) E. 
e 



D„ = 



2e 



d s (n) Ap, 



(22) 
(23) 



where D a has a linear dependence on the average elec- 
tric field through a proportionality constant d a which 
depends on the doping n. D s is instead the contribution 
to the interlayer polarization coming from the width of 
the transferred charge. Therefore we write it in a form 
consistent with the other interlayer term in Eq.(18), i.e., 



de 2 /(2eo)Ap, with a proportionality constant d s which 
depends on the doping n. 

In Fig jlO] we show d a and d s as a function of doping 
n. Since d s is almost independent of the doping and 
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FIG. 11: (Color online) a(n) calculated with DFT-GGA (x 
crosses) from Eq.(12l using the DFT-GGA calculated rj(n) 



(squares) and from Eq. ( 25 1 using the DFT-GGA calculated 
rj(n) (stars). 



FIG. 12: (Color online) The DFT-GGA calculated a as a 
function of the doping n for an electronic temperature of 300 
(crosses) and 30 K (stars). 



d a has a variation in the order of 5%, we replace them 
with their average values calculated on the doping range 
considered, d„ = l-09 A, and J s =0.80 A. Using only this 



approximation and Eq.(13l, we have 



U=(3 



-e(d - d a )E & 



2e 



(d-d s )Ap 



(24) 



We notice that the simplified electrostatic model de- 
scribed in Sec IIB i.e., Eq.(|6|, used in TB calculations, 
is equivalent to consider, in Eq.(24l, (3 — 1, d a = 

, p^{z) = 5{\z\ - 



{i.e. p^a\z) — 0), and d s 
d/2) Ap/2 sign(z)]. 







Considering Eq.(24l, and the definition of a [as in 
Eq.Q] and rj [as in Eq.Q] we obtain another relation 
between a and rj as follows: 



a(n) — 



°/3 (d-d a )/d 



1 - 77(71) a baro /3 (d-d s )/d 



(25) 



Equation ( 25 1 gives the approximate relation between 



U, the average electric field, the screening charge ob- 
tained considering the intralaycr polarization, and con- 
sidering the width of the transferred charge between lay- 
ers. This equation substitutes Eq.( 12 ) which comes from 



the simplified electrostatic model described in Sec jll B 

In Fig 11 we show the DFT-calculated a(n), a(n) ob- 
tained from Eq.(12l using the DFT-calculated i](n), and 
from the electrostatic model of Eq.(25l using the DFT- 



calculated rj{n). One can see that the simplified elec- 
trostatic model is not able to describe the DFT results. 



Instead, cv(n) obtained from the model of Eq.(25) is able 
to correctly reproduce the DFT calculations. 



Effect of the electronic temperature on a as a 
function of doping n 



In Sec. |III A| we have shown that in the bilayer the 
electronic screening to the external electric field is cru- 
cial for a correct evaluation of the band-gap. At low 
doping level the screening is expected to depend on the 
broadening parameter, and in this section we investigate 
the effect of the electronic temperature on the screening 
and on a. 



In Fig 12 we show the DFT-calculated a as a function 
of n for an electronic temperature of 300 and 30 K. The 
variation in screening with the broadening parameter de- 
pends on the doping n. Since the doping levels which are 
interesting for applications of the bilayer as active de- 
vice in nanoelectronics are small values around the zero 
doping, we focus on this doping range. 



In Fig 13 we show the DFT-calculated a as a function 
of the electronic temperature T for electron doping values 
between and 5.72 x 10 12 cm~ 2 . In this range of doping, 
we can see that the difference between a at 300 and 30 K 
is largest for zero doping. In particular, for zero doping 
the band-gap at 30 K results to be about 10% smaller 
than at 300 K. 



D. Nonmonotonic behavior of a as a function of 
doping n 

As shown in Figs(5]and|6] both DFT-calculated a and 
•q have a nonmonotonic behavior as a function of the 
doping n. a and rj represent the linear response of U 
to the external average electric field E av and the linear 
response of Ap to the band-gap U, respectively. Up to 
now, we calculated a and rj for finite values of E av and U. 
In this section we show that perturbation theory (PT) 
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FIG. 13: (Color online) The DFT-GGA calculated a as a 
function of the electronic temperature T, for small values 
of electron doping. The values of doping n are in units of 

ml2 -2 



FIG. 14: (Color online) Schematic representation of allowed 
(dashed arrows) and not allowed (continuous arrows) contri- 
butions from band-to-band transitions in Eq.( 28 1, for electron 
doping. 



explains the origin of this nonmonotonic behavior as a 
function of doping n. 

For the numerical evaluation of the expressions ob- 
tained from PT, we use the band structure calculated 
with the TB model. Indeed, even if TB results for a 
and 77 differ from the DFT ones, TB is able to catch the 
nonmonotonic trend of these quantities as a function of 
the doping n. Moreover, we limit our PT calculations 
to 77(77). Indeed, since the relation between 77 and a is 
monotonic [see Eq.(25)], the nonmonotonic behavior of 



77 is able to explain also the nonmonotonic behavior of a 
as a function of n. 

In order to calculate 77(71) with PT, we consider 

which is the unperturbed TB Hamiltonian. is a 

4x4 matrix which depends on the wave vector k and is 
written on the basis of 2p z orbitals centered on the four 
atoms of the unit-cell, ordered as A, B, A', B' (A and B 
are the two carbon atoms on layer 1, A' and B' are the 
two carbon atoms on layer 2, and in the Bernal stacking 
configuration A and A' are vertically superposed). In 
presence of a band splitting U (see Fig|2|, the Hamilto- 
nian i/k can be written as 



(o) 



U 



(26) 



where 



Ap = 



1 














1 














-1 

















(27) 



Using first-order PT, we obtain the following expres- 
sionfor77= (^): 



1 



(0) 
k 



,(0) 



JO) _ (0) 



(0) 



1 (k 



\<^\Ap\4Z> 



> 



(28) 



where 



> (7=1,2,3,4) are the unperturbed eigen- 



' ' °' ) with eigenvalues 



states of if k 

/(e|k — £ f' 1 ) i s the occupation of state i 



is the Fermi level, 



e' F ' j is tne occupation ot state 7, and Nk is the 
number of k points used in the BZ integration. 

In particular, in Eq.(28) there are six contributions 



>, 



obtained by mixing the four unperturbed states \ip^ 
with 7=1,2,3, and 4; and we label as r/ij the contribu- 
tion to 77 given by states i and j. Within the present 
TB model 7713 and 7724 are exactly zeroJ^S Contribution 
771,2 vanishes for ep > because the two states are both 
occupied. Therefore, for 6f > 0, the important contribu- 
tions to 77 derive from 771,4, 772,3, and 773,4, as schematically 
shown in Fig [14] 



In Fig 15 we show 77 as a function of the electron dop- 
ing, obtained from Eq.(28), for an electronic temperature 
of 300 K. Different contributions from 771,4, 772,3, and 773 4 
are also plotted. For comparison, we report 77(71) calcu- 
lated with the nonperturbative TB model. Contribution 
771,4 as a function of doping is constant when the Fermi 
level is lower than the bottom of band 4, and its absolute 
value starts to decrease when band 4 becomes occupied 
due to lower availability of empty states. Contribution 
i]2,3 is minimum for zero doping, and its absolute values 
decreases, as a function of the electron doping, for the 
same reason. Contribution from 77^4 is lower than con- 
tribution from 772,3 due to the energy difference in the 
denominator of Eq.(28l, which is higher for 771,4. Finally, 
contribution 773,4 vanishes for zero doping, and its abso- 
lute value increases with increasing electron doping, due 
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FIG. 15: (Color online) 77 as a function of electron doping, 
from Eq.(28l (linear response), and different contributions 
from 771,4, 772,3, and 773,4- Squares are the nonperturbative 
TB results. 



to larger number of possible transitions between occupied 
and unoccupied states. 



The results in Fig 15 show that the nonmonotonic be- 
havior of 77(71) is determined by the sum of the three con- 
tributions 771,4, ??2,3) and 773,4, which gives a maximum at 
the electron doping values n rs 6 x 



E. GW correction 

Recently it has been shown by ARPES measure- 
ments that the electronic band structure of graphene and 
graphite is not well reproduced by LDA and GGAP3 In 
particular, LDA and GGA underestimate the slope of the 
bands since these approximations do not include long- 
range electron-correlation effects. Such effects can fully 
be taken into account within the GW approach [where 
the self-energy is computed from the product of the elec- 
tron Green's function (G) and the screened Coulomb in- 
teraction (W)], which is considered to be the most ac- 
curate first principles approach for the electronic band 
structure.^ The GW band structures for graphite and 
graphene are indeed in very good agreement with the 
ARPES measurements.^ In absence of an external elec- 
tric field, the DFT-calculated bands of bilayer graphene 
need to be scaled in order to reproduce the GVF-correct 
bands as 



C GW 



\ ,DFT 



(29) 



where A = 1.18 is the scaling factor, as obtained from 
RefJ^S. Such scaling factor can change the screening 
properties of the bilayer, and in this section we include 
it in our theoretical results. 

If we focus on the quantity 77, we can use the per- 
turbative expression in Eq.(28). In such expression, we 



TABLE I: TB-GW 7 parameters obtained by fitting the bilayer 
DFT bands with a TB model, along all the FKM line, and by 
rescaling the parameters with A = 1.18. 7! is the i-nearest- 
neighbors hopping parameters. All values are in eV. 





7j} 

















TB-GW -3.4013 0.3292 -0.2411 0.1226 0.0898 0.3963 0.1671 0.3301 



correct the DFT eigenvalues using Eq.(|29|, and we 
can neglect the GW correction to the matrix elements 



A/9 IV'lk^ ^ smcc it is commonly found that the 



DFT error in the wave functions is usually negligible with 
respect to the error on the eigenvalues. Within this ap- 
proximation, it is easy to show that 



DFTV^F ? 

1 A ' A 



(30) 



where T is the temperature. The computed i] GW is 
shown in Fig (6] 

a GW can be computed from i] GW using our model in 
Eq.(25). In order to minimize the error from our model, 



we write 



a 



(«) 



3 /3(d-d„)/d 



1 - r/ GW (n) a halc f3(d-d s )/d 



A, (31) 



where 



DFT 



„barc 



(3(d-d a )/d 



1 -77 DFT a barc ^(d- d s )/d' 



(32) 



gives an estimate of the error in Eq.(|25|. The computed 
a GW (n) is shown in Figj5] and for low doping levels it is 
around 10% higher than the DFT value. 



F. Full band structure of gated bilayer graphene 

In this section we give a practical instruction to obtain 
the full band structure of bilayer graphene for a doping n 
and for an average electric field E av . In order to do that, 
we fit our DFT bands along all the TKM line in the BZ 
in absence of the external electric field, using a TB model 
with five nearest neighbors in-plane hopping parameters 
^\\ ' ' ' ^\\ ' ) an< ^ ^hree out-of-plane hopping param- 
eters (7^ ,J± B iJ± B )• In the Bernal stacking config- 
uration of bilayer graphene, A and A' represent the ver- 
tically superposed atoms. These hopping parameters do 
not change when an external average electric field is ap- 
plied. This is shown in Fig |16[ where we compare the 
direct DFT results with the TB calculations, with the 
fixed hopping parameters and the U value from the DFT 
calculations. 

Since we consider the GW one as the most precise re- 
sult, in Table |T] we report the TB-GW hopping parame- 
ters obtained by fitting the DFT bands without electric 
field and by rescaling them with the GW scaling factor 
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FIG. 16: (Color online) DFT-GGA calculated bands around FIG. 17: (Color online) Comparison between experimental 
the K point in the BZ (DFT) for n = and with U = 0.14 eV results for U from RefP (squares), TB, and GW results, 
and from the TB model (TB fit), with five nearest-neighbors 
in-plane hopping parameters and three out-of-plane hopping 
parameters. 

G. Comparison with experimental results 



TABLE II: Values of fitting parameters of Eq.( 33 1. All values 
are in 10 -12 cm 2 meV. 



Ai 


0.896 


Si 


-26.888 


7i 


21.756 


A 2 


3.905 


B 2 


1.623 


72 


21.946 


A 3 


-1.654 


B 3 


-0.092 


73 


5.534 


C 


5.848 



A = 1.18. Moreover, in order to avoid the numerical eval- 
uation of U for a given n and E av , we give a fit of our 
calculated a GW (n): 



a GW (n) 



3 

E 

i=l 



(n-B z ) 2 
if 



C. 



(33) 



where the values of the fitting parameters are listed in 
Table [IT] In Fig(5] we show the results of the fit with 
the black continuous line. From expression (33 1, we can 



obtain the value of the gap U as a function of the doping 
n and of the external average electric field E av , 



U(n,E av ) = a GW {n){ ni -n 2 ), 



(34) 



where (nx — n 2 ) = -E av /(|e|/(2e )). n, n\, and n 2 are in 
units of 10 12 cm~ 2 . 



In this section we compare our DFT and TB results 
with the direct measurement (ARPES) of the band-gap 
on epitaxially growth bilayer graphene in RefP. In this 
work Ohta et al. 3 performed an experiment where bi- 
layer graphene is synthesized on silicon carbide (SiC) sub- 
strate. The SiC acts as a fixed bottom gate, and a charge 
ri2 flows from the substrate to the bilayer. Further elec- 
tron doping is induced with the deposition of potassium 
atoms on the other side of the bilayer, and this chemi- 
cal doping acts as a top gate. Varying the concentration 
of potassium, the asymmetry between the two layers of 
graphene is modified, and a band-gap is opened accord- 
ingly. Using angle-resolved photoemission spectroscopy 
Ohta et aiP directly measured the band structure, and 
fitting it with a TB model, they obtained a curve of the 
gap as a function of the doping charge in the bilayer. 

To compare with their experimental res ults, wecalcu- 
late the gap using a GW (n) from Eqs.(33l and (34 1 and 
keep n 2 (bottom gate) fixed at 11.9 x 10 12 cm~ . This 
value of n.2 derives from the fact that in RefP, for a to- 
tal doping of n = 23.8 x 10 12 cm~ 2 , no gap is observed, 
meaning that ri\ — n-i — n/2. Since in the experiment 
the bottom gate is not varied, we also keep it fixed to 
this value, and we only vary n\ = n — 

we compare our results, obtained with a 



In Fig 



17 



and with a TB , with the experimental data from RefP 
We first notice that the nonlinearity is not due to the 
saturation of the gap with E& v ; it is instead due to the 
dependence of a on the doping n (at high doping a de- 
creases with n). Moreover, both GW and TB results are 
in good agreement with experiments. This is due to the 
fact that the experiment is carried out at high doping 
levels, where the difference between the GW and TB a's 
is less important with respect to low doping levels (see 
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Fig§. 

In the case of exfoliated bilayer graphene, direct ex- 
perimental measurements of the band structure and of 
the gap with ARPES are still unavailable. Alternatively, 
indirect information on the band structure can be ob- 
tained by infrared reflectivity studies. Recently, Kuz- 
menko et aZP^ reported an experimental work on in- 
frared spectra of exfoliated and gated bilayer graphene 
as a function of doping. In this work the authors found a 
strong gate- voltage dependence of their spectral features, 
which are related to interband transitions. A compari- 
son of the experimental infrared spectra with the one ob- 
tained from TB calculations suggests that the TB predic- 
tion of gate-induced band-gap is overestimated.^ How- 
ever, a quantitative analysis of the band-gap as a function 
of doping and external field is not given. 

Finally, by measuring the cyclotron mass as a function 
of doping in bilayer graphene one can check the pres- 
ence of a finite band-gap. These measurements do not 
provide a direct estimate of the band-gap; however, they 
give important informations on the hole-electron asym- 
metry and on the deformation of the band structure in 
the presence of an external electric field. In particular, 
in Refs.^ ancP the authors measured the cyclotron mass 
on exfoliated bilayer graphene. The bottom gate is re- 
alized with an oxidized silicon substrate, which allows a 
variation in bottom gate electron density ri2 during the 
experiment. The top gate is provided by chemical dop- 
ing, by deposition of NH3 molecules, which provides a 
top gate electron density n 1; which is fixed during the 
experiment. 

To compare the experimental results of Refs.^ ancP 
with our band structures, we calculate the cyclotron mass 
m r as 



m c (n) 



r ( dA{E) \ 
27 \ dE J 



(35) 



E=E f (n) 



where A is the k-space area enclosed by the orbit with 
energy E and Ef is the Fermi level. The derivative in 
Eq.(35l is obtained by finite differentiation with respect 



to E. For the GW calculations we use the band structure 
calculated as described in Sec. IIIIFl 

In Fig. [l8]we compare the experim enta l data on the cy- 
clotron mass from Refplwith our [Fig 18 'a)] GW calcula- 
tions and with [Fig 18 'b)] TB calculations^ for different 
values of top gate electron density n\. In RefP the au- 
thors estimated an initial doping no on bilayer graphene, 
at zero bottom gate, of about 1.8 xl0 12 cm~ 2 . In prin- 
ciple, such initial doping could come both from the de- 
posited NH3 molecules {i.e., from the top gate) and from 
a charge transfer from the SiC"2 substrate (i.e., from the 
bottom gate). Thus an exact estimation of the top gate 
electron density n\ is not possible, and we calculate the 
cyclotron mass for values of ri\ between 1.8 xl0 12 cm -2 
and 0. Our results show that for both GW and TB cal- 
culations, the cyclotron mass behavior as a function of 
doping depends on the value of n\. In particular, for the 
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FIG. 18: (Color online) Cyclotron mass, with respect to the 
free electron mass m e , as a function of doping n: comparison 
of the experimental results from RefP with our (a) GW cal- 
culations and with our (b) TB calculations for different values 
of ni (the values of m are in units of 10 12 cm -2 ). 



GW calculations the best agreement with the experimen- 
tal results is obtained for m = 0.45 xl0 12 cm -2 . Finally, 
we note that our GW calculations give better results than 
the TB ones. In particular, contrary to the TB results, 
they arc able to reproduce the hole-electron asymmetry. 



IV. CONCLUSIONS 

We present a detailed ab initio DFT investigation of 
the band-gap opening and screening effects of gated bi- 
layer graphene. First, we analyze the response of the 
band-gap to the external average electric field at fixed 
doping. We show that this response is linear for differ- 
ent electron and hole doping values and for large electric 
field values. We then find that the linear response of the 
gap to the electric field has a nonmonotonic behavior as 
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Supercell 



FIG. 19: Planar average of the ionic, Hartree, and monopole 
potentials, multiplied by the electron charge. This fig- 
ure corresponds to a doping charge on the bilayer n = 
19 x 10 12 cm" 2 . The positions of the first and second lay- 
ers of the bilayer and of the monopole in the supercell are 
indicated. 



a function of doping and for low doping values it depends 
on the temperature. 

We also perform TB calculations for the band-gap 
opening. At low doping values, which are the interesting 
ones for electronic applications, we find that the DFT- 
calculated gap is roughly half of the TB one. Since the 
band-gap strongly depends on the screening effects, we 
perform a detailed analysis of the charge distribution in 
the bilayer in presence of the external electric field. We 
show that the electronic screening is characterized by in- 
terlayer and intralayer polarizations. The latter one, not 
included in TB calculations, gives an important contri- 
bution to the band-gap opening. 

On the basis of this analysis, we propose a model which 
significantly improves the description of the electronic 
properties of bilayer graphene in the presence of an ex- 
ternal electric field, and finally we provide a practical 
scheme to obtain the full band structure of gated bilayer 
graphene for arbitrary values of the doping and of the 
external electric field. 
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APPENDIX A: DIPOLE AND MONOPOLE 
POTENTIAL 

Standard plane-wave ab initio codes work with period- 
ically repeated super-cells. When doping the sample with 



a total charge neA (A is the area of the section of the 
periodic cell parallel to the graphene plane) a compen- 
sating uniform background charge (with opposite sign) 
is added in order to have a neutral system and, thus, a 
periodic electrostatic potential. 

In this work, we use a different approach, and we add 
a "monopole" , that is a uniformly charged plane equidis- 
tant from the two graphene layers, with total charge 
—neA. This is done by adding in real space a periodic 
potential energy given by 



?7,e 2 z 2 
V mon (z) = -— (-z+ — ), 



(Al) 



where z = z — z mon , z mon is the z coordinate of the 
monopole plane, and z £ [0',L], where L is the length 
of the supercell along z. While the linear term of V" m0 n 
is the potential associated with the monopole plane, the 
quadratic term cancels the potential associated with uni- 
form background. V mon can be added to the electrostatic 
potential acting on the Kohn-Sham electronic states with 
a straightforward implementation. The resulting system 
is, as a whole, neutral and periodic. 
In Fig. 



19 we show the planar average of the ionic, 



Hartree, and monopole potentials, multiplied by the elec- 
tron charge. The position of the first and second layers of 
the bilayer in the supercell is indicated, together with the 
monopole position. The distance between the monopole 
and the bilayer is 6.93 A. This figure corresponds to a 
doping charge on the bilayer n = 19 x 10 12 cm" 2 and to 
an experimental setup where the bottom and top gates 
are equal, and no gap opening is expected. 

In order to have different bottom and top gates, we 
add to the monopole a sawlike potential, called dipole 
potential,^ generated by two planes of opposite charge, 
as implemented in standard distributions of the PWSCF 
code^M The dipole is centered around the monopole, 
and the distance bet wee n the dipole planes is kept fixed 
to 0.17 A. 



In Fig. 



20 



we show the planar average of 
the ionic, Hartree, monopole, and dipole potentials mul- 
tiplied by the electron charge for a doping charge n = 
19 x 10 12 cm" 2 . In the case shown in the figure, the 
dipole potential is chosen to create a flat potential and 
zero electric field on layer 1 of the bilayer. This configu- 
ration corresponds to the case where only a bottom gate 
acts on the bilayer. By changing the sign to the dipole 
potential, we can obtain the opposite configuration, with 
a flat potential and zero electric field on layer 2 of the 
bilayer. 

The electric fields E\ and E 2 are calculated from the 
planar average of the ionic, Hartree, monopole, and 
dipole potential energy V\{z) and V2(z) on side 1 and 
side 2 of the bilayer, respectively, as 



E 2 = 



1\ dV^z) 
e / dz ' 
1\ dV 2 (z) 
el dz 



(A2) 
(A3) 
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Layer 1 Layer 2 Monopole Layer 1 




In order to deal with uniform E\ and E 2 electric fields, 
these derivatives are calculated in the linear part of V\(z) 
and V 2 (z) (see Figj20j. 

Varying independently the dipolc potential and the to- 
tal charge on the sample and monopole, one can explore 
all the situations with different doping n on the bilayer 
and different E av = (E 1 + E 2 )/2. 



Supercell 



FIG. 20: Planar average of the ionic, Hartree, monopole, and 
dipole potentials multiplied by the electron charge. The po- 
sitions of the first and second layers of the bilayer, of the 
monopole, and of the dipole in the supercell are indicated. 
This figure corresponds to a doping charge on the bilayer n 
= 19 x 10 12 cm -2 . The dipole potential is such that layer 1 
of the bilayer does not feel any external electric field. 
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